// stud.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"

using std::cout;
using std::endl;

using boost::format;

using std::string;

using std::copy;

#define DIM(a) (sizeof(a) / sizeof(a[0]))
struct V2 {
	float x;
	float y;

	V2(float x = 0.0f, float y = 0.0f) : x(x), y(y) {
	}
};

struct V3 {
	float x;
	float y;
	float z;

	V3(float x = 0.0f, float y = 0.0f, float z = 0.0f) : x(x), y(y), z(z) {
	}
};

using std::vector;

typedef vector<V2> V2Vec;
typedef vector<V3> V3Vec;

void create_stud_section(V2Vec& vec, float size, float flange,
	float thickness, float lip){
	vec.clear();

	V2 section[] = {
		V2(0.0f, 0.0f), V2(0.0f, size), V2(flange, size),
		V2(flange, size - lip), V2(flange - thickness, size - lip),
		V2(flange - thickness, size - thickness),
		V2(thickness, size - thickness), V2(thickness, thickness),
		V2(flange - thickness, thickness),
		V2(flange - thickness, lip), V2(flange, lip), V2(flange, 0.0f)
	};

	vec.insert(vec.begin(), section, section + sizeof(section) / sizeof(V2));
}

template<typename T>
int push_back(vector<T>& vec, const T v){
	int rv = vec.size();
	vec.push_back(v);
	return rv;
}

const int QSEG = 8;

void create_circle(V2Vec& circle, float radius){
	V2Vec qne, qse, qsw, qnw;

	for(int i = QSEG - 1; i >= 0; --i){
		double theta = i * M_PI / 2.0 / QSEG;
		double stheta = sin(theta);
		double ctheta = cos(theta);

		float h = (float)(radius * ctheta);
		float v = (float)(radius * stheta);

		qne.push_back(V2(h, v));
		qse.push_back(V2(v, -h));
		qsw.push_back(V2(-h, -v));
		qnw.push_back(V2(-v, h));
	}

	circle.clear();
	circle.insert(circle.end(), qnw.begin(), qnw.end());
	circle.insert(circle.end(), qne.begin(), qne.end());
	circle.insert(circle.end(), qse.begin(), qse.end());
	circle.insert(circle.end(), qsw.begin(), qsw.end());

}

string polystring(const vector<int>& poly){
	string& s = str(format("%d ") % poly.size());
	for(vector<int>::const_iterator i = poly.begin(), end = poly.end();
		i != end; ++i)
		s = s + str(format("%d ") % *i);

	return s;
}

void report(const V3Vec& points,
		const vector<int>& fnear, const vector<int>& ffar,
		const vector<int>& hhard, const vector<int>& hsoft);

void init_stud(V3Vec& points, vector<int>& fnear, vector<int>& ffar,
	vector<int>& hhard, vector<int>& hsoft){

	points.clear();
	fnear.clear();
	ffar.clear();
	hhard.clear();
	hsoft.clear();

	const float size = 6.0f;
	const float flange = 1.625f;
	const float thickness = 0.033f;
	const float length = 119.934f;
	const float lip = 0.5f;

	// hole dimensions
	const float h_width = 1.0f;
	const float h_height = 3.0f;
	const float h_radius = h_width / 2.0f;
	const float h_position = 36.0f;

	V2Vec section;
	create_stud_section(section, size, flange, thickness, lip);

	for(V2Vec::iterator itr = section.begin(); itr != section.end(); ++itr)
		fnear.push_back(push_back(points, V3(itr->x, itr->y, 0.0f)));

	for(V2Vec::iterator itr = section.begin(); itr != section.end(); ++itr)
		ffar.push_back(push_back(points, V3(itr->x, itr->y, length)));

	V2Vec circle;
	create_circle(circle, h_radius);

	assert(circle.size() % 2 == 0);

	V2Vec uarc, larc;
	int n = circle.size();
	uarc.insert(uarc.begin(), circle.begin(), circle.begin() + n / 2);
	larc.insert(larc.begin(), circle.begin() + n / 2, circle.end());

	float offset = (h_height - 2 * h_radius) / 2.0f;
	for(int i = 0, m = n / 2; i < m; i++){
		V2& u = uarc[i];
		u.y += offset;

		V2& l = larc[i];
		l.y -= offset;
	}

	V2Vec& hole = circle;
	hole.clear();

	for(V2Vec::iterator i = uarc.begin(), end = uarc.end(); i != end; ++i)
		hole.push_back(*i);

	for(V2Vec::iterator i = larc.begin(), end = larc.end(); i != end; ++i)
		hole.push_back(*i);

	float holey = size / 2.0f;
	for(V2Vec::iterator i = hole.begin(), end = hole.end(); i != end; ++i)
		hhard.push_back(push_back(points,
			V3(0.0f, holey - i->x, h_position + i->y)));

	for(V2Vec::iterator i = hole.begin(), end = hole.end(); i != end; ++i)
		hsoft.push_back(push_back(points,
			V3(thickness, holey - i->x, h_position + i->y)));

}

void report(const V3Vec& points,
		const vector<int>& fnear, const vector<int>& ffar,
		const vector<int>& hhard, const vector<int>& hsoft){

	cout << format("%d 3 0 0") % points.size() << endl;
	for(int i = 0, n = points.size(); i < n; i++){
		const V3& v = points[i];
		cout << format("%d %f %f %f") % i % v.x % v.y % v.z << endl;
	}

	cout << format("%d 0") % (2 + ffar.size() + hsoft.size()) << endl;

	cout << "1 0" << endl;
	cout << polystring(fnear) << endl;

	cout << "1 0" << endl;
	cout << polystring(ffar) << endl;

	for(int i = 0, n = fnear.size(), lasti = n - 1; i < n; ++i){
		int nexti;
		if(i == lasti)
			nexti = 0;
		else
			nexti = i + 1;

		if(i == 0){
			cout << "2 1" << endl;
			int face[4] = { fnear[i], fnear[nexti], ffar[nexti], ffar[i] };
			cout << polystring(vector<int>(face, face + DIM(face))) << endl;
			cout << polystring(hhard) << endl;
			V3 p0 = points[hhard[0]], pmid = points[hhard[hhard.size() / 2]];
			V3 vhole((p0.x + pmid.x) * 0.5f, (p0.y + pmid.y) * 0.5f,
				(p0.z + pmid.z) * 0.5f);
			cout << format("1 %f %f %f") % vhole.x % vhole.y % vhole.z << endl;
		}
		else if(i == 6){
			cout << "2 1" << endl;
			int face[4] = { fnear[i], fnear[nexti], ffar[nexti], ffar[i] };
			cout << polystring(vector<int>(face, face + DIM(face))) << endl;
			cout << polystring(hsoft) << endl;
			V3 p0 = points[hsoft[0]], pmid = points[hsoft[hsoft.size() / 2]];
			V3 vhole((p0.x + pmid.x) * 0.5f, (p0.y + pmid.y) * 0.5f,
				(p0.z + pmid.z) * 0.5f);
			cout << format("1 %f %f %f") % vhole.x % vhole.y % vhole.z << endl;
		}
		else{
			cout << "1 0" << endl;
			int face[4] = { fnear[i], fnear[nexti], ffar[nexti], ffar[i] };
			cout << polystring(vector<int>(face, face + DIM(face))) << endl;
		}
	}

	for(int i = 0, n = hhard.size(), lasti = n - 1; i < n; i++){
		int nexti;
		if(i == lasti)
			nexti = 0;
		else
			nexti = i + 1;

		cout << "1 0" << endl;
		int face[4] = { hhard[i], hhard[nexti], hsoft[nexti], hsoft[i] };
		cout << polystring(vector<int>(face, face + DIM(face))) << endl;

	}

	cout << "0" << endl;
	cout << "0" << endl;
	
}

void init_polygon(tetgenio::polygon& poly, const vector<int>& verts){
	poly.numberofvertices = verts.size();
	poly.vertexlist = new int[poly.numberofvertices];
	copy(verts.begin(), verts.end(), poly.vertexlist);
}

void init_polygon(tetgenio::polygon& poly, const int *begin, const int *end){

	poly.numberofvertices = end - begin;
	poly.vertexlist = new int[poly.numberofvertices];
	copy(begin, end, poly.vertexlist);
}

void init_tetgen(tetgenio& tgio, const V3Vec& points,
	const vector<int>& fnear, const vector<int>& ffar,
	const vector<int>& hhard, const vector<int>& hsoft){

	tgio.firstnumber = 0;
	tgio.numberofpoints = points.size();
	tgio.pointlist = new REAL[tgio.numberofpoints * 3];
	REAL *p = tgio.pointlist;
	for(int i = 0, n = points.size(); i < n; i++){
		const V3& v = points[i];
		*p = v.x;
		*(p + 1) = v.y;
		*(p + 2) = v.z;
		p += 3;
	}

	tgio.numberoffacets = 2 + fnear.size() + hhard.size();
	tgio.facetlist = new tetgenio::facet[tgio.numberoffacets];
//	tgio.facetmarkerlist = new int[tgio.numberoffacets];
//	for(int i = 0, n = tgio.numberoffacets; i < n; ++i)
//		tgio.facetmarkerlist[i] = 42;

	tetgenio::facet *pfacet = tgio.facetlist;

	pfacet->numberofpolygons = 1;
	pfacet->polygonlist = new tetgenio::polygon[pfacet->numberofpolygons];
	init_polygon(pfacet->polygonlist[0], fnear);
	pfacet->numberofholes = 0;
	pfacet->holelist = 0;
	pfacet++;

	pfacet->numberofpolygons = 1;
	pfacet->polygonlist = new tetgenio::polygon[pfacet->numberofpolygons];
	init_polygon(pfacet->polygonlist[0], ffar);
	pfacet->numberofholes = 0;
	pfacet->holelist = 0;
	pfacet++;

	for(int i = 0, n = fnear.size(), lasti = n - 1; i < n; i++){
		int nexti;
		if(i == lasti)
			nexti = 0;
		else
			nexti = i + 1;

		int face[4] = { fnear[i], fnear[nexti], ffar[nexti], ffar[i] };

		if(i == 0){
			pfacet->numberofpolygons = 2;
			pfacet->polygonlist =
				new tetgenio::polygon[pfacet->numberofpolygons];
			init_polygon(pfacet->polygonlist[0], &face[0], &face[4]);
			init_polygon(pfacet->polygonlist[1], hhard);

			V3 p0 = points[hhard[0]], pmid = points[hhard[hhard.size() / 2]];
			V3 vhole((p0.x + pmid.x) * 0.5f, (p0.y + pmid.y) * 0.5f,
				(p0.z + pmid.z) * 0.5f);
			pfacet->numberofholes = 1;
			pfacet->holelist = new REAL[3];
			REAL *phole = pfacet->holelist;
			phole[0] = vhole.x; phole[1] = vhole.y; phole[2] = vhole.z;
		}
		else if(i == 6){
			pfacet->numberofpolygons = 2;
			pfacet->polygonlist =
				new tetgenio::polygon[pfacet->numberofpolygons];
			init_polygon(pfacet->polygonlist[0], &face[0], &face[4]);
			init_polygon(pfacet->polygonlist[1], hsoft);

			V3 p0 = points[hsoft[0]], pmid = points[hsoft[hsoft.size() / 2]];
			V3 vhole((p0.x + pmid.x) * 0.5f, (p0.y + pmid.y) * 0.5f,
				(p0.z + pmid.z) * 0.5f);
			pfacet->numberofholes = 1;
			pfacet->holelist = new REAL[3];
			REAL *phole = pfacet->holelist;
			phole[0] = vhole.x; phole[1] = vhole.y; phole[2] = vhole.z;
		}
		else {
			pfacet->numberofpolygons = 1;
			pfacet->polygonlist =
				new tetgenio::polygon[pfacet->numberofpolygons];
			init_polygon(pfacet->polygonlist[0], &face[0], &face[4]);
			pfacet->numberofholes = 0;
			pfacet->holelist = 0;
		}

		pfacet++;
	}

	for(int i = 0, n = hhard.size(), lasti = n - 1; i < n; i++){
		int nexti;
		if(i == lasti)
			nexti = 0;
		else
			nexti = i + 1;

		int face[4] = { hhard[i], hhard[nexti], hsoft[nexti], hsoft[i] };

		pfacet->numberofpolygons = 1;
		pfacet->polygonlist =
			new tetgenio::polygon[pfacet->numberofpolygons];
		init_polygon(pfacet->polygonlist[0], &face[0], &face[4]);
		pfacet->numberofholes = 0;
		pfacet->holelist = 0;

		pfacet++;
	}

	tgio.pointmarkerlist = new int[tgio.numberofpoints];
	for(int i = 0, n = tgio.numberofpoints; i < n; i++){
		tgio.pointmarkerlist[i] = 32;
	}

	assert(tgio.numberoffacets == (pfacet - &tgio.facetlist[0]));
}

int _tmain(int argc, _TCHAR* argv[])
{

	//V2Vec circle;
	//create_circle(circle, 1.0f);

	//for(V2Vec::iterator i = circle.begin(); i != circle.end(); ++i){
	//	cout << format("%f, %f") % i->x % i->y << std::endl;
	//}

	V3Vec points;
	vector<int> fnear, ffar, hhard, hsoft;
	init_stud(points, fnear, ffar, hhard, hsoft);
	
	//report(points, fnear, ffar, hhard, hsoft);

	tetgenio tg, tgout;
	init_tetgen(tg, points, fnear, ffar, hhard, hsoft);
	tetrahedralize("pYze", &tg, &tgout);

	cout << "==================" << endl;

	int *ptri = tgout.trifacelist;
	int nmypoints = 0;
	for(int i = 0, n = tgout.numberofpoints; i < n; ++i)
		if(tgout.pointmarkerlist[i] == 32)
			nmypoints++;
		else
			cout << "unknown " << tgout.pointmarkerlist[i] << endl;

	cout << "nmypoints = " << nmypoints << endl;
	//for(int i = 0, n = tgout.numberoftrifaces; i < n; ++i){
	//	cout << str(format("%d %d %d") % ptri[0] % ptri[1] % ptri[2]) << endl;
	//	ptri += 3;
	//}


	return 0;
}
